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This report summarizes the effort expended during the first six 

months on the study of supersonic turbulent boundary layers over two- 

dimensional protuberances. The approach taken is to extend an earlier 

work (Ref. 1) to the turbulent regime, using the nximerical finite- 

difference alternating direction implicit (ADI) method, fm obvious 

departure from the previous case (Ref, 1) is forced by the need to 

model mathematically the turbulence. The turbulence is represented 

here by the eddy viscosity approach as used in Refs. (2) and (3). 

The turbulent boundary layer structure as well as an interest in thick 

boundary layers and much larger protuberance heights than in the 

laminar case lead to new difficulties. The problems encountered in 

the course of the present investigation and the means to remove them 

will be discussed later in this report. 

The governing equations pertinent to the probl^ are 

Continuity 

V + F + 2^F, = 


0 
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Momentum 

3A 2 

(£eF^)^ “ + (l+a/2) (s+ ^ (g-P ) - 2ePP, = 0 

n n n o u ^ 

Energy 

A 

(- 1 ^)^ - Vg^ [^(e-e/Pr) FP^]^ - 25Pg^ = 0 

with 

F(5,0) = V(e,0) = 0 , g(4,0) = Hyn^ 

F(5,“) = g(5/“) = 1 . 

Ths momentum equation is cast into the present form to facilitate the 
use of the ADI algorithm. The notation is the same as in Refs. (1) 
and (3) , Hence, F = normalized streamwise component of velocity, 
g = normalized total enthalpy, V = transformed normal velocity function 
= total displacement body, I = molecular viscosity parameter, a = hn 

^ A 

inviscid parameter, and e and e are eddy viscosity parameters defined 
for fully turbulent flow as e = 1 + e/y, e = 1 + = — s/yr where e = 
eddy viscosity. The Prandtl numbers, Pr and Pr^ are taken to be 
constant. 

The study commenced with programming the above equations using 

the Cebeci-Smith two-layer eddy viscosity model, applied successfully 

* 

to attached turbulent boundary layers earlie.r (Refs. (2), (3)), A 


* 


Initially, this choice is adequate since the immediate aim is to 
develop an efficient nume’^ical algorithm; comparison with experi- 
mental data will be delayed for later. 
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base case was identified with flow conditions corresponding to the 

NASA test data (Ref. (4)}. The free stream Mach number and static 

temperature are M^ = 2,5 and = 252°R respectively. The Reynolds 

* 

number based on free stream conditions and a reference length L = 

15.25 cm is Re^ = 1.647x10^, The sine-wave protuberance profile is 
given by y = ^ {1 - cos (x - 1 ) t where h is the height and 

w is the width of the protuberance. In all the base case studies 
w = 0.24 (i.e. w* = 0.24 X 15.25 cm = 3.66 cm; star quantities are 
dimensional). For a single wave placed on a flat plate far downstreaun 
from the leading edge, at s = 25.85, and h = 0.02 a calculation was 
performed with T^T^ 0.81, Pr = 0.72 and Pr^ = 0,9. An attached 
solution was obtained (Fig. 1). In this case the boundary layer 
displacement thickness 6 is about three times the protuberance height 
(but only about 1/3 as thick as the NASA Langley tunnel wall boundary 
layer). The calculation commenced at s^ = 1.0 with a laminar boundary 
layer, passing through transition and becoming fully turbulent at 
s = 1.4. After sweeping only once in %, direction up to station s^ = 
25.5 the ADI algorithm is employed between stations s^ and s^ = 26.5, 

The protuberance is placed between x = 25.85 and 26.09. The streaun- 
wise step size. As was taken 0.02 in this calculation so that there 
are at least 12 grid points over the protuberance surface. In 6.7 
minutes of computer time 22 time sweeps were performed and the skin 
friction is varying (between the last two time steps at a given station) 
in the third or higher decimal place. This was considered a converged 
solution. It is worthwhile noticing that in the laminar case (Ref. 1) 
the same number of time steps required about 15-20 minutes computer 
time. The present efficiency was accomplished by linearization of the 
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coefficients of the corresponding difference equation around the 
previous station value, eliminating thus the need to update these 
coefficients in subsequent iterations. 

Increasing the protuberance height at otherwise identical flow 
conditions as given above, a separation region appeared ahead of the 
protuberance. But the values of skin friction (the pareuneter we 
monitor for testing the quality of the solution) did not converge in 
time. On the contrary, wide oscillations occured in time indicating 
some numerical instability, and eventually the calculation terminated 
abruptly. Examination of the boundary layer velocity profiles showed 
slight oscillations (of order 10 ^ around the value P = 1) in the 
boundary layer edge region. At this point it was not clear whether 
this was due to a stronger interaction (higher protuberance) , the 
reverse flow region, an improper eddy model or some other effects. 

To isolate the source of trouble it was decided to simplify the problem 
by taking Pr *= Pr^ = 1, thus eliminating the need to solve the energy 
equation. (This corresponds to adiabatic wall condition with g = 
constant = 1) ; the rationale being that since the artificial time de- 
pendent term appears only in the momentum equation, the source of the 
time-like oscillations maybe due to finite-difference representation 
of this equation. Because the difficulty appeared for the separated 
cases, the eddy viscosity model became suspect as well. 'Freezing' the 
eddy viscosity distribution (this idea was borrowed from Ref. 5) near 
the forward plate-protuberance junction, and the use of Alber's model 
(Ref. 6) for the separated region did not resolve the difficulty. 

It was brought to our attention (Ref, 7) that upwind differencing of 
the term in the momentiam equation maybe required to satisfy the 
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convergence criteria of the numerical scheme (see also Ref, 8), In 
the present case the momentum equation, after linearization, may be 
written as 

F + a,F + ct-F + a, + ay,F_ + Oe V = 0 , 

nnln234c5 ' 

where a^, ... Og are the linearized coefficients. In the boundary 

layer near the wall the diffusion term F is predominant over the 

nn 

convection- like term o.F and a central difference scheme for P 

In n 

is appropriate. From a so-called diagonal dominance test of the model 
equation F + aF = 0 it is found that with central differencing the 
criteria |aAn| 2 must be met. Because of additional terms in. the 
momentum equation we have found (by numerical experimentation) that 
the use of the criteria |aAni 1 eliminated the oscillations in the 
profile near the boundary layer edge. Hence, F^ is central differenced 
when |aAn| ^ 1 and upwind differencing is used when jaAnj > 1. Davis 
(Ref. 7) also pointed out that the use of forward differencing of F^ 
term in the continuity equation is more appropriate because the upstream 
propagation becomes important for the strongly interacting boundary 
layers. Having incorporated these ideas into our algorithm difficulties 
mentioned above have been removed; the boundary layer profile approaches 
smoothly the edge values (F = 1) and the time-like oscillations also 
disappeared. For the base flow conditions, with the forward junction 
of the plate-protuberance placed at s^ = 3.35 and protuberance height 
of 0.024 results are shown in Figures 2 and 3. The dashed line in 
Figure 2 represents the variation of wall shear with time steps at a 
fixed station (s = 3.33) located just downstream of the plate-protu- 
berance junction. This typically unstable solution, obtained before 
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the above mentioned changes were introduced ^ terminates abruptly after 
ten time steps. After the modifications in the numerical algor ittua 
were made an apparently convergent solution was obtained (full line. 

Fig. 2) . The variation in appears only in the third or higher 

decimal place after the first fifteen time steps. The spatial dis- 
tribution of at the last time step (33rd sweep) is plotted in Figure 
3. A small separation bubble appears near the forward plate-protuberance 
junction and the maximvun shear is about double that of the flat plate 

value. The boundary layer displacement thickness ahead of the protu- 

* * 

berance is relatively thin, 6 0.010 (or 6 = 5L = 0.1525 cm) and 

6/h = 0.4. 

A major interest of this research project is in comparisons of 

the numerical predictions with experimental data of very thick separated 

turbulent boundary layers produced on the tunnel wall (Ref. 4). The 

separation characteristics over the test plate placed in the wind tunnel 

wall are strongly dependent on the approaching turbulent boundary layer 

profile. For this reason calculations were performed for the ^ jundary 

layer as it develops along the wall of the UPWT Langley Tunnel. The 

following test section free-stream conditions were taken; = 2.535, 

Reycm = 1.08x10^, T^ = 567°R (315^K) , T^ = 252°R, T^T^ = 0.81, 

* 

Pr = 0.72, Pr- = 0.9. A reference length L = 15.25 cm (= 0.5 ft). 

^ p*u*L* 

Correspondingly, the reference Reynolds number is Re = 9 P — = 

oo 

1.08x10^x15,25 = 1.647x10^; and the nondimen sional protuberance width 
w of 0.24 corresponds to w of 3.66 cm. The boundary layer calculation 
was carried out as a non-interacting 2-D laminar-transitional-turbulent 
boundary layer developing from ahead of the nozzle throat under a 
favorable pressure gradient. In the supersonic region downstream of 
the throat the pressure distribution was calculated from the sidewall 
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Mach number distribution, using isentropic relations. The Mach 
number distribution was obtained from a characteristic net (Ref. 9). 
From these data a cubic representation for M was assumed of the fom 


M 


“t ■ 


K, 


(1 - f--) 


(1 - K 


s 
2 s. 


where test section Mach number at s » s^. The above polynomial 

representation satisfies the condition dM/ds » 0 at s = s^. At the 

location s of the first characteristic near the throat the Mach number 
o 

was estimated to be M « 1.11. The measured distance along the wall 

o 

from s to sm is 39.6 ('x* 20 ft). At a location s_ 8.86 units downstream 

Ct X (5 

of s the Mach number is M. * 1.84. Letting s =1, s„=9.86 and s- 
ot p a p X 

= 40.6. Using these values, the constants and K 2 can be determined. 
For K^, K 2 > 0 and K 2 ^ 1 (which is the case here) the polynomial rep- 
resentation given above yields monotomically increasing Mach number. 

The subsonic-transonic section of the tunnel was also assumed to be 
represented by a cubic polynomial 


M - M. + (1 - 1^) (1 - Cj I-) , 

with the following properties: 

At s = s. = 0.3, M = M. and dM/ds = 0. 

1 X 

At s . s M = M = 1.11 and ^ ^ 

s 

Three different values for M^ were chosen: 0.01, 0.03 and 0,05, 

It turned out that the boundary layer development downstream of the 
throat v'as not sensitive to these initial values. Taking M^ = 0.03 
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a boundary layer calculation was performed for 0.3 ^ s ^ 97.4 with 
nine c<^plete boundary layer profiles punched on IBM cards between 
locations 40.6 and 72.6. (At s » 40.6 the Mach number becomes con- 
stant on the side wall; at s ^ 50.6 the straight section begins). 

These profiles are presently being used as initial data for separation 
studies of thick turbulent boundary layers. The table belcv^ gives the 
calculated boundary layer displacement thickness distribution along 
the constant Idach number section of the UPWT Langley tunnel wall at 
ten stations 


s 40.6 

44.6 

48.6 

52.6 

56.6 

60.6 

64.6 

68.6 

•72,6 

76.6 

6* (cm) 1.69 

1.84 

1.98 

2.12 

2.26 

2.40 

2.54 

2.67 

2.80 

2.93 


Using the profile produced at station s = 72.60 and placing a 
single sine-wave protuberance on the flat with its junction at = 

73.25 calculations were initiated to obtain separation characteristics 
of thick turbulent boundary layers. The calculation starts at s * 72.60 
by sweeping once to station 72.96 and then the 7U)I algorithm was employed 
between this station and s = 73.96 using 51 grid points in the s direct- 
ion (As = 0.02) with 105 point variable mesh in n direction. Figure 
4 shows results obtained for the base case flow conditions with h « 

0.08, w = 0,24, Pr = 0.72, Pr^ « 0.90 and T^T^ = 0.81. Note that here 
the energy equation was solved simultaneously with the continuity and 
momentum equations. Figure 4 shows the variation of the surface skin 
friction and pressure levels predicted by the present method. It is 
seen that ahead of the protuberance the skin friction drops and a 
small separation region develops on the front face of the wave. The 
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skin friction peaks just aft of the wave top with a value three times 
the undisturbed flat plate value. The pressure distribution is some- 
what different and shows a monctomic increase up to a peak near the 
top of the wave followed by a mild recovery back to the flat plate 

value. 

While the above results demonstrc'te the capability of the present 
method to handle flow fields typical of present interest, there re- 
main several refinements to be made before the technique can be con- 
sidered riperational. 

The first of these involves the size of the finite difference 
mesh applied to the problem. Using the present mesh from s = 73 to 
74 it is noticed that the surface pressure distribution for all points 
in the mesh ahead of an aft of the protuberance is influenced by the 
protuberance. This effect is manifest by the depressed pressure levels 
(p/p^ < 1) near s = 73 and the monotomically increasing levels of 
ahead of protuberance. A test was conducted with the ends of the finite 
difference mesh moved away from the protuberance. This lead to a re- 
lief of this difficulty. Apparently a larger number of grid points 
ahead of the protuberance will have to be used. 

The second, area of interest needing further study is the conver- 
gence rate of the iterative scheme. The present ADI numerical algorithm 
is conditionally stable and the convergence rate sensitive to the time 
step chosen. For the solutions depicted in Figure 4, convergence was 
achieved in about 40 iterations (approximately 10 minutes) and this 
computer time can be further reduced by either reducing the number of 
grid points (presently set at 5000) or parametrically assessing the 
optimum value of At. Through these two efforts it would be anticipated 
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that the computer effort could be reduced to the 2 or 3 minute level. 
Efforts in this direction will be made during the second half of this 
study. 

Another point of interest involves the sensitivity of the pre- 
sent solution method to the initialization procedure. At present 
the solution method is found to show some influence of initial con- 
ditions (time B 0) on the final steady state values, mainly in the 
separated r ion and the source of this anomaly is presently not known. 
This problem could be due to either a programming error (a possibility 
now being investigated) or a conceptual error. It is possible that 
the influence of the confined finite difference mesh discussed above 
is causing this problem, but it is not yet clear how this would occur. 
Another possibility is that the consistency error of the ADI scheme 
(Ref. 10) is producing this difference and this possiblity will also 
be considered. 

Following this an eddy viscosity model which is more appropriate 
for separated boundary layers will be incorporated into the algorithm 
before comparisons with experimental data for single and multiple 
protuberances will be made. 
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